setwd('/data1/bsi/bioinf_int/s106381.borawork/data_HTHGU/maqc/affy_plus2_taqman/whole_rma/1')
library(affy)
library(frma)
data<-ReadAffy()
library(simpleaffy)
library(frma)
library(hgu133plus2frmavecs)
library(MASS)
library(preprocessCore)



frma <- function(object, background="rma", normalize="quantile", summarize="robust_weighted_average", input.vecs=list(normVec=NULL, probeVec=NULL, probeVarBetween=NULL, probeVarWithin=NULL, probesetSD=NULL), output.param=NULL, verbose=FALSE){

  if(!is(object, "AffyBatch")) stop(paste("argument is", class(object), "frma requires AffyBatch"))

  if(!background %in% c("none", "rma")) stop("background must be either none or rma")
  if(!normalize %in% c("none", "quantile")) stop("normalize must be either none or quantile")
  if(!summarize %in% c("median_polish", "average", "median", "weighted_average", "robust_weighted_average", "batch")) stop("summarize must be one of: median_polish, average, median, weighted_average, robust_weighted_average, batch")

  if(summarize=="batch" & (any(input.vecs$probeVarBetween==0) | any(input.vecs$probeVarWithin==0))) stop("If summarize method is batch then probeVarBetween and probeVarWithin must be greater than zero for all probes.")
  
  cdfname <- cleancdfname(cdfName(object))
  platform <- gsub("cdf","",cdfname)

  if(summarize == "median_polish") output <- frmaMedPol(object, background, normalize, input.vecs, verbose)
  if(summarize %in% c("average", "median", "weighted_average", "robust_weighted_average")) output <- frmaRobReg(object, background, normalize, summarize, input.vecs, output.param, verbose)
  if(summarize == "batch") output <- frmaBatch(object, background, normalize, input.vecs, output.param, verbose)

  if("stderr" %in% output.param) stderr <- output$stderr else stderr <- matrix(nrow=0, ncol=0)
  if("weights" %in% output.param) w <- output$weights else w <- matrix(nrow=0, ncol=0)
  if("residuals" %in% output.param) r <- output$residuals else r <- matrix(nrow=0, ncol=0)

  R.model <- PLM.designmatrix3(object)
  
  new("PLMset",
      chip.coefs=output$exprs,
      weights=list("PM.weights"=w, "MM.weights"=matrix(nrow=0, ncol=0)),
      se.chip.coefs=stderr,
      residuals=list("PM.resid"=r, "MM.resid"=matrix(nrow=0, ncol=0)),
      cdfName = object@cdfName,
      phenoData = phenoData(object),
      annotation = annotation(object),
      experimentData = experimentData(object),
      nrow= object@nrow,
      ncol= object@ncol,
      narrays=length(object),
      model.description = list("R.model"=R.model))
}
frmaBatch <- function(object, background, normalize, input.vecs, output.param, verbose){

  cdfname <- cleancdfname(cdfName(object))
  platform <- gsub("cdf","",cdfname)
  
  if(background == "rma"){
    if(verbose) message("Background Correcting ...\n")
    object <- bg.correct.rma(object)
  }

  if(is.null(input.vecs$normVec) | is.null(input.vecs$probeVec) | is.null(input.vecs$probeVarWithin) | is.null(input.vecs$probeVarBetween)){
    pkg <- paste(platform, "frmavecs", sep="")
    require(pkg, character.only=TRUE, quietly=TRUE) || stop(paste(pkg, "package must be installed first"))
    data(list=eval(pkg))

    if(is.null(input.vecs$normVec)) input.vecs$normVec <- get(pkg)$normVec
    if(is.null(input.vecs$probeVec)) input.vecs$probeVec <- get(pkg)$probeVec
    if(is.null(input.vecs$probeVarWithin)) input.vecs$probeVarWithin <- get(pkg)$probeVarWithin
    if(is.null(input.vecs$probeVarBetween)) input.vecs$probeVarBetween <- get(pkg)$probeVarBetween
  }

  if(normalize == "quantile"){
    if(verbose) message("Normalizing ...\n")
    pm(object) <- normalize.quantiles.use.target(pm(object), input.vecs$normVec)
  }

  if(verbose) message("Summarizing ...\n")
  pns <- probeNames(object)
  pms <- log2(pm(object))

  N <- 1:dim(pms)[1]
  S <- split(N, pns)

  fit <- lapply(1:length(S), function(i) {
    s <- S[[i]]
    batchFit(pms[s,], input.vecs$probeVec[s], input.vecs$probeVarWithin[s], input.vecs$probeVarBetween[s])
  })
  
  exprs <- matrix(unlist(lapply(fit, function(x) x$exprs)), ncol=ncol(pms), byrow=TRUE)
  rownames(exprs) <- names(fit)
  colnames(exprs) <- colnames(pms)

  if("weights" %in% output.param){
    weights <- matrix(unlist(lapply(fit, function(x) x$w)), ncol=ncol(pms), byrow=TRUE)
    rownames(weights) <- pns
    colnames(weights) <- sampleNames(object)
  } else weights <- NULL

  if("stderr" %in% output.param){
    stderr <- unlist(lapply(fit, function(x) x$se))
    names(stderr) <- names(fit)
  } else stderr <- NULL

  if("residuals" %in% output.param){
    residuals <- apply(exprs,2,function(x) rep(x, table(pns)))
    residuals <- (pms-input.vecs$probeVec) - residuals
    rownames(residuals) <- pns
    colnames(residuals) <- sampleNames(object)
  } else residuals <- NULL
  
  colnames(exprs) <- sampleNames(object)
  
  return(list(exprs=exprs, stderr=stderr, weights=weights, residuals=residuals))
}

batchFit <- function(x1, x2, x3, x4){
  y.new <- x1 - x2
  yy.new <- as.vector(t(y.new))
  X <- matrix(rep(diag(ncol(y.new)),nrow(y.new)), nrow=length(yy.new), byrow=TRUE)
  Z <- diag(nrow(y.new))
  Z <- Z[rep(seq(nrow(Z)), each=ncol(y.new)),]
  G <- diag(x4)
  R <- diag(rep(x3, each=ncol(y.new)))
  V <- Z%*%G%*%t(Z) + R
  Vinv <- solve(V)
  e <- eigen(Vinv, symmetric=TRUE)
  C <- (e$vectors)%*%diag(sqrt(e$values))%*%t(e$vectors)
  yy.trans <- C%*%yy.new
  x.trans <- C%*%X

  fit.batch <- rlm(yy.trans ~ -1 + x.trans, maxit=100)
  b <- fit.batch$coef
  names(b) <- NULL
  u <- G%*%t(Z)%*%Vinv%*%(yy.new-X%*%b)
  offset.u <- mean(u)
  exprs <- b + offset.u
  w <- fit.batch$w
  se <- fit.batch$s
  return(list(exprs=exprs, w=w, se=se))
}
frmaRobReg <- function(object, background, normalize, summarize, input.vecs, output.param, verbose){

  cdfname <- cleancdfname(cdfName(object))
  platform <- gsub("cdf","",cdfname)
  
  if(background == "rma"){
    if(verbose) message("Background Correcting ...\n")
    object <- bg.correct.rma(object)
  }
  
  if(is.null(input.vecs$normVec) | is.null(input.vecs$probeVec) | is.null(input.vecs$probeVarWithin) | is.null(input.vecs$probeVarBetween) | (summarize=="robust_weighted_average" & is.null(input.vecs$probesetSD))){
    pkg <- paste(platform, "frmavecs", sep="")
    require(pkg, character.only=TRUE, quiet=TRUE) || stop(paste(pkg, "package must be installed first"))
    data(list=eval(pkg))

    if(is.null(input.vecs$normVec)) input.vecs$normVec <- get(pkg)$normVec
    if(is.null(input.vecs$probeVec)) input.vecs$probeVec <- get(pkg)$probeVec
    if(is.null(input.vecs$probeVarWithin)) input.vecs$probeVarWithin <- get(pkg)$probeVarWithin
    if(is.null(input.vecs$probeVarBetween)) input.vecs$probeVarBetween <- get(pkg)$probeVarBetween
    if(is.null(input.vecs$probesetSD) & summarize=="robust_weighted_average") input.vecs$probesetSD <- get(pkg)$probesetSD
  }

  if(normalize == "quantile"){
    if(verbose) message("Normalizing ...\n")
    pm(object) <- normalize.quantiles.use.target(pm(object), input.vecs$normVec)
  }

  if(verbose) message("Summarizing ...\n")
  pns <- probeNames(object)
  pms <- log2(pm(object))
  
  if(summarize == "average"){  
    exprs <- subColSummarizeAvg((pms-input.vecs$probeVec), pns)
    weights <- NULL
    stderr <- NULL
  }

  if(summarize == "median"){  
    exprs <- subColSummarizeMedian((pms-input.vecs$probeVec), pns)
    weights <- NULL
    stderr <- NULL
  }
  
  if(summarize == "weighted_average"){  
    w <- 1/(input.vecs$probeVarWithin + input.vecs$probeVarBetween)
    if(any(input.vecs$probeVarWithin==0) | any(input.vecs$probeVarBetween==0)) warning("Either probeVarWithin or probeVarBetween is 0 for some probes -- setting corresponding weights to 1")
    w[w==Inf] <- 1
    exprs <- subColSummarizeAvg((pms-input.vecs$probeVec)*w, pns)
    W <- as.vector(rowsum(w, pns, reorder=FALSE))
    exprs <- (exprs/W)*as.vector(rowsum(rep(1,length(pns)),pns,reorder=FALSE))
    weights <- NULL
    stderr <- NULL
  }

  if(summarize == "robust_weighted_average"){
    w <- 1/(input.vecs$probeVarWithin + input.vecs$probeVarBetween)
    if(any(input.vecs$probeVarWithin==0) | any(input.vecs$probeVarBetween==0)) warning("Either probeVarWithin or probeVarBetween is 0 for some probes -- setting corresponding weights to 1")
    w[w==Inf] <- 1
    N <- 1:dim(pms)[1]
    S <- split(N, pns)
    fit <- lapply(1:length(S), function(i) {
	s <- S[[i]]
	rwaFit2(pms[s,, drop=FALSE], w[s], input.vecs$probeVec[s], input.vecs$probesetSD[i])
    })
    names(fit) <- unique(pns)
    exprs <- matrix(unlist(lapply(fit, function(x) x$Estimates)), ncol=ncol(pms), byrow=TRUE)
    rownames(exprs) <- names(fit)
    colnames(exprs) <- colnames(pms)
    if("weights" %in% output.param){
      weights <- matrix(unlist(lapply(fit, function(x) x$Weights)), ncol=ncol(pms), byrow=TRUE)
      rownames(weights) <- pns
      colnames(weights) <- sampleNames(object)
    } else weights <- NULL
    if("stderr" %in% output.param){
      stderr <- matrix(unlist(lapply(fit, function(x) x$StdErrors)), ncol=ncol(pms), byrow=TRUE)
      rownames(stderr) <- names(fit)
      colnames(stderr) <- sampleNames(object)
    } else stderr <- NULL
  }

  if("residuals" %in% output.param){
    residuals <- apply(exprs,2,function(x) rep(x, table(pns)))
    residuals <- (pms-input.vecs$probeVec) - residuals
    rownames(residuals) <- pns
    colnames(residuals) <- sampleNames(object)
  } else residuals <- NULL
  
  colnames(exprs) <- sampleNames(object)
  
  return(list(exprs=exprs, stderr=stderr, weights=weights, residuals=residuals))
}

rwaFit2 <- function(x1, x2, x3, x4){
  ncols <- ncol(x1)
  w.tmp <- x2/max(x2)
  w.tmp <- matrix(rep(w.tmp, ncols), ncol=ncols)
  pe.tmp <- x3
  pe.tmp[1] <- pe.tmp[1]-sum(pe.tmp)
  rcModelWPLM(y=x1, w=w.tmp, row.effects=pe.tmp, input.scale=x4)
}

barcode <- function(object, platform=NULL, mu=NULL, tau=NULL, cutoff=6.5, output="binary"){

  if(!class(object) %in% c("matrix", "ExpressionSet", "PLMset") & !is.vector(object)) stop("Object must be one of: vector, matrix, ExpressionSet, or PLMset.")

  if(!output %in% c("binary", "z-score", "p-value", "lod")) stop("Output must be one of: binary, z-score, p-value, lod.") 

  if(is.vector(object)) object <- as.matrix(object)
  
  if(is.matrix(object)){
    if(is.null(platform)) stop("If object is of class matrix or vector, platform cannot be NULL.")
    if(!platform %in% c("GPL96", "GPL570", "GPL1261")) stop("Platform must be one of: GPL96, GPL570, GPL1261")
  }

  if(class(object)=="ExpressionSet"){
    if(annotation(object) %in% c("hgu133a", "hgu133plus2", "mouse4302")){
      if(cleancdfname(annotation(object)) == "hgu133acdf") platform <- "GPL96"
      if(cleancdfname(annotation(object)) == "hgu133plus2cdf") platform <- "GPL570"
      if(cleancdfname(annotation(object)) == "mouse4302cdf") platform <- "GPL1261"
    } else stop("Microarray platform given by object annotation not recognized -- please supply platform type.")
    object <- as.matrix(exprs(object))
  }

  if(class(object)=="PLMset"){
    if(annotation(object) %in% c("hgu133a", "hgu133plus2", "mouse4302")){
      if(cleancdfname(annotation(object)) == "hgu133acdf") platform <- "GPL96"
      if(cleancdfname(annotation(object)) == "hgu133plus2cdf") platform <- "GPL570"
      if(cleancdfname(annotation(object)) == "mouse4302cdf") platform <- "GPL1261"
    } else stop("Microarray platform given by object annotation not recognized -- please supply platform type.")
    object <- as.matrix(coefs(object))
  }

  i.remove <- grep("AFFX", rownames(object))
  if(length(i.remove)>0) object <- as.matrix(object[-i.remove,])

  if(platform == "GPL96" & nrow(object)!=22215) stop("Object does not have the correct dimensions for platform GPL96.")
  if(platform == "GPL570" & nrow(object)!=54613) stop("Object does not have the correct dimensions for platform GPL570.")
  if(platform == "GPL1261" & nrow(object)!=45037) stop("Object does not have the correct dimensions for platform GPL1261.")
  
  if(is.null(mu) | is.null(tau)){
    if(platform=="GPL96") pkg <- "hgu133abarcodevecs"
    if(platform=="GPL570") pkg <- "hgu133plus2barcodevecs"
    if(platform=="GPL1261") pkg <- "mouse4302barcodevecs"
    require(pkg, character.only=TRUE, quiet=TRUE) || stop(paste(pkg, "package must be installed first"))
    data(list=paste("bcparams-", platform, sep=""))
    if(is.null(mu)) mu <- bcparams[,2]
    if(is.null(tau)) tau <- sqrt(bcparams[,3])
  }
  
  e <- object
  num <- round(nrow(e)/10)
  
  if(output %in% c("p-value", "lod", "binary")){
    pval <- pnorm(e, mean=mu, sd=tau, lower.tail=FALSE)
    if(output == "p-value"){
      colnames(pval) <- colnames(e)
      rownames(pval) <- rownames(e)
      return(pval)
    } else{
      lod <- -log10(pval)
      if(output == "lod"){
        colnames(lod) <- colnames(e)
        rownames(lod) <- rownames(e)
        return(lod)
      } else{
        bc <- matrix(as.integer(lod > cutoff), ncol=ncol(e))
        colnames(bc) <- colnames(e)
        rownames(bc) <- rownames(e)
        return(bc)
      }
    }
  }

  if(output == "z-score"){
    z <- sweep(sweep(e, 1, mu), 1, tau, FUN="/")
    colnames(z) <- colnames(e)
    rownames(z) <- rownames(e)
    return(z)
  }
}














fr<-frma(data,summarize="batch")
write.exprs(fr,"normalized_new_frma.txt")
